﻿// intersection_mesh.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include <iostream>
#include <string>
#include "mesh.hpp"
#include "GlobalMainVars.h"


extern void readVTK(Mesh& mesh, const std::string& filename);

int mesh()
{
    // std::string filename = "./OUTPUT.vtk";
	cout << "Reading mesh hexahedron..." << endl;
	std::string filename = "./hexahedron.vtk";
    HexMesh hexMesh;
    hexMesh.readVTK(filename);
    hexMesh.initialize();
    //filename = "./Tetrahedralization.vtk";
	cout << "Reading mesh Tetrahedralization..." << endl;
	filename = "./Tetrahedralization.vtk";
    //TetMesh tetMesh;
    tetMesh.readVTK(filename);
	tetMesh.initialize();
	cout << "Finished reading and initializing the mesh!" << endl;
	cout << "Performing interpolation..." << endl;
	size_t count = 0;
	cout << "Number of elements in tetMesh: " << tetMesh.elems.size() << endl;
#pragma omp parallel for num_threads(8)
	for (int i = 0;i < tetMesh.elems.size();++i) {
		auto& elem = tetMesh.elems[i];
		if (elem->rockType != 1) continue;
		auto tetElem = dynamic_cast<TetElement*>(elem.get());
		Point p = { tetElem->center };
		int idx = hexMesh.find_containing_cell(p);
		try {
			/*if (idx == -1) {
				cerr << "Point not in mesh!" << endl;
				cerr << "count = " << count << endl;
				cerr << "Point: " << p << endl;
			}*/
			// now point p is inside hexMesh.elems[idx], use the value for interpolation
			elem->data = hexMesh.elems.at(idx)->data;
			++count;
			if (count % 1000 == 0) {
				cout << "Interpolated " << count << " elements" << endl;
				//cout << "idx = " << idx << endl;
			}
		}
		catch (std::out_of_range) {
			cerr << "Point not in mesh!" << endl;
			cerr << "count = " << count << endl;
			cerr << "Point: " << p << endl;
			//break;
		}
	}
	return 0;
}